# Load Necessary Packages
library(glmnet)
library(foreign)
library(estimatr)
library(dplyr)
library(tidyr)
library(bcf)
library(dbarts)
library(texreg)
library(Hmisc)
library(grid)
library(gridExtra)
library(interplot)
library(margins)
library(stringr)
library(factoextra)
library(ggplot2)
library(stargazer)
library(TOSTER)
library(interflex)
library(effectsize)

## Load Data
# Load Data
setwd("/users/josephphillips/Dropbox/COVID-19/Data/")
us <- read.dta("2021 01 - Granite State Poll - Weighted (Dartmouth) v12.dta")

## Load glmnet function
run_LASSO <- function(data, yvar, pre){
  
  set.seed(2334)
  
  var_quo <- enquo(yvar)
  var_quo1 <- enquo(pre)
  
  ## select x variables 
  dat <- data %>%
    dplyr::select(college, agegroup, male, married, frequent_church, pid3, ideo7, lives_highincidence, nonwhite, health_trust)
  
  pre <- data %>%
    dplyr::select(!!var_quo1) %>% pull() # Save as a vector instead of a data frame ot use model matrix
  
  ## select y-variable
  out1 <- data %>% 
    dplyr::select(!!var_quo) %>% pull() # save as a vector instead of a data-frame to use model.matrix()
  
  merged1 <- cbind(out1, dat, pre) %>% na.omit() ## merge and delete NA
  x1 <- model.matrix(out1 ~ ., data = merged1) ## make model matrix
  y1 <- merged1$out1 ## select y-variable from NA-deleted dataset
  mod1 <- cv.glmnet(x1, y1, alpha = 1, family = "gaussian") ## run lasso model
  coef <- coef(mod1) ## get coef  
  
  list(merged1, x1, y1, mod1, coef)
  
}

# Clean variables
us$treatment <- factor(us$SPLIT10,levels=c("Control","Injunction","Contention"))
us$dv_missing <- ifelse(us$Q602=="One Dose" | us$Q602=="Two Doses",NA,7-as.numeric(us$COVID19_Vaccine))
us$dv_verylikely <- ifelse(us$Q602=="One Dose" | us$Q602=="Two Doses",6,7-as.numeric(us$COVID19_Vaccine))
us$dv_mostlikely <- ifelse(us$Q602=="One Dose" | us$Q602=="Two Doses",7,7-as.numeric(us$COVID19_Vaccine))
us$college <- ifelse(us$DEMOS_education=="College graduate" | us$DEMOS_education=="Postgraduate work",1,0)
us$agegroup <- us$DEMOS_age
us$male <- ifelse(us$DEMOS_gender=="Men",1,0)
us$married <- ifelse(us$DEMOS_marital=="Married",1,0)
us$frequent_church <- ifelse(us$DEMOS_religious_services=="Once a week or more" | us$DEMOS_religious_services=="Once or twice a month",1,0)
us$pid3 <- as.numeric(us$DEMOS_partyid)
us$ideo7 <- as.numeric(us$DEMOS_ideology)
us$lives_highincidence <- ifelse(us$DEMOS_county=="Hillsborough",1,0)
us$nonwhite <- ifelse(as.numeric(us$DEMOS_race_caucasian_white)==2,0,1)
us$health_trust <- ifelse(us$COVID7=="Don't trust",0,ifelse(us$COVID7=="Unsure",1,2))
us$lagdv_missing <- ifelse(us$Q602=="One Dose" | us$Q602=="Two Doses",NA,ifelse(us$Q438=="Almost certainly not",0,ifelse(us$Q438=="Probably not",1,ifelse(us$Q438=="Probably",2,ifelse(us$Q438=="Almost certainly",3,NA)))))
us$lagdv_verylikely <- ifelse(us$Q602=="One Dose" | us$Q602=="Two Doses",3,ifelse(us$Q438=="Almost certainly not",0,ifelse(us$Q438=="Probably not",1,ifelse(us$Q438=="Probably",2,ifelse(us$Q438=="Almost certainly",3,NA)))))
us$lagdv_mostlikely <- ifelse(us$Q602=="One Dose" | us$Q602=="Two Doses",4,ifelse(us$Q438=="Almost certainly not",0,ifelse(us$Q438=="Probably not",1,ifelse(us$Q438=="Probably",2,ifelse(us$Q438=="Almost certainly",3,NA)))))

# Lasso Procedure
df <- us
lasso_missing <- run_LASSO(df,dv_missing,lagdv_missing) ## pre
lasso_verylikely <- run_LASSO(df,dv_verylikely,lagdv_verylikely) ## pre
lasso_mostlikely <- run_LASSO(df,dv_mostlikely,lagdv_mostlikely) ## pre

# Balance, Table S1
table(us$treatment,us$college)
chisq.test(us$treatment,us$college)
table(us$treatment,us$agegroup)
us$age1834 <- ifelse(us$agegroup=="18 to 34",1,0)
us$age3549 <- ifelse(us$agegroup=="35 to 49",1,0)
us$age5064 <- ifelse(us$agegroup=="50 to 64",1,0)
us$age65 <- ifelse(us$agegroup=="65 and older",1,0)
chisq.test(us$treatment,us$age1834)
chisq.test(us$treatment,us$age3549)
chisq.test(us$treatment,us$age5064)
chisq.test(us$treatment,us$age65)
table(us$treatment,us$male)
chisq.test(us$treatment,us$male)
table(us$treatment,us$married)
chisq.test(us$treatment,us$married)
table(us$treatment,us$frequent_church)
chisq.test(us$treatment,us$frequent_church)
table(us$treatment,us$pid3)
us$dem <- ifelse(us$pid3==1,1,0)
us$rep <- ifelse(us$pid3==3,1,0)
chisq.test(us$treatment,us$dem)
chisq.test(us$treatment,us$rep)
mean(subset(us,treatment=="Control")$ideo7,na.rm=T)
mean(subset(us,treatment=="Injunction")$ideo7,na.rm=T)
mean(subset(us,treatment=="Contention")$ideo7,na.rm=T)
summary(aov(ideo7~treatment,data=us))
table(us$treatment,us$lives_highincidence)
chisq.test(us$treatment,us$lives_highincidence)
table(us$treatment,us$nonwhite)
chisq.test(us$treatment,us$nonwhite)
mean(subset(us,treatment=="Control")$health_trust,na.rm=T)
mean(subset(us,treatment=="Injunction")$health_trust,na.rm=T)
mean(subset(us,treatment=="Contention")$health_trust,na.rm=T)
summary(aov(health_trust~treatment,data=us))
mean(subset(us,treatment=="Control")$lagdv_verylikely,na.rm=T)
mean(subset(us,treatment=="Injunction")$lagdv_verylikely,na.rm=T)
mean(subset(us,treatment=="Contention")$lagdv_verylikely,na.rm=T)
summary(aov(lagdv_verylikely~treatment,data=us))
mean(subset(us,treatment=="Control")$lagdv_mostlikely,na.rm=T)
mean(subset(us,treatment=="Injunction")$lagdv_mostlikely,na.rm=T)
mean(subset(us,treatment=="Contention")$lagdv_mostlikely,na.rm=T)
summary(aov(lagdv_mostlikely~treatment,data=us))
mean(subset(us,treatment=="Control")$lagdv_missing,na.rm=T)
mean(subset(us,treatment=="Injunction")$lagdv_missing,na.rm=T)
mean(subset(us,treatment=="Contention")$lagdv_missing,na.rm=T)
summary(aov(lagdv_missing~treatment,data=us))

# Table S6/Figure 1
Reg.verylikely.us.robust <- lm_robust(dv_verylikely~treatment+lagdv_verylikely,data=us)
us1 <- effectsize(Reg.verylikely.us.robust)
Reg.mostlikely.us.robust <- lm_robust(dv_mostlikely~treatment+lagdv_mostlikely,data=us)
Reg.missing.us.robust <- lm_robust(dv_missing~treatment+lagdv_missing,data=us)

# Table S14
Reg.verylikely.us.robust.nocov <- lm_robust(dv_verylikely~treatment,data=us)
Reg.mostlikely.us.robust.nocov <- lm_robust(dv_mostlikely~treatment,data=us)
Reg.missing.us.robust.nocov <- lm_robust(dv_missing~treatment,data=us)